Code
#|label: setup
#|echo: false
knitr::opts_chunk$set(
message = FALSE, warning = FALSE, collapse = TRUE
)#|label: setup
#|echo: false
knitr::opts_chunk$set(
message = FALSE, warning = FALSE, collapse = TRUE
)This dataset is come from Koval and his colleagues’ (2013) study in which they use a novel paradigm and anlytical approach to study the dynamic process of depression. In past study, depression is associated with higher average levels of negative affect (NA) and lower average levels of positive affect (PA). However, few studies untangle the relationship between the pattern of affective fluctuations and the dynamics of depression. Therefore, in their study, they collected 99 (the final number) subjects’ depression symptom by the Center for Epidemiologic Studies Depression Scaleand (CES-D) and daily experiences of affect ratings 10 times/day for 7 days using the experience sampling method (ESM). Affect ratings were made on two items measuring PA (= mean(happy, relaxed)) and four items measuring NA (=mean(sad, depressed, anxious, angry)), with the scale from 1 (not at all) to 100 (very much).
#|label: load-packages-and-rawdata
library(tidyverse)
library(lubridate)
rawdata <- read_tsv("data_1beep_no1st beep_annette.tsv")
rmarkdown::paged_table(rawdata)#|label: sorting-data
data <- rawdata %>%
mutate(Pos = PA,
Neg = `NA`,
Date_Time = ymd_hms(str_glue("{year}-{month}-{day} {hour}:{min}:{sec}"), tz = "CET"),
Date = as_date(Date_Time),
Time = hms::as_hms(Date_Time),
WDay = wday(Date, label = TRUE),
Subject = factor(cumsum(PpID != lag(PpID, default = 0))),
.keep = "none") %>%
group_by(Subject) %>%
mutate(Day = factor(cumsum(Date != lag(Date, default = origin)))) %>%
group_by(Subject, Date) %>%
mutate(Moment = factor(1:n())) %>%
ungroup() %>%
pivot_longer(cols = c(Pos, Neg), names_to = "Affect", values_to = "Score")
rmarkdown::paged_table(data)#|label: line-plot-of-score-in-diff-day-time
data %>%
filter(Subject %in% c(6, 17, 44, 60)) %>%
drop_na() %>%
ggplot(aes(x = as_datetime(paste(today(), Time)),
y = Score,
group = WDay, color = WDay)) +
geom_point() +
geom_line() +
scale_x_datetime(date_breaks = "2 hour", date_labels = "%H:%M",
limits = c(as.POSIXct(c("9:30", "23:30"), format = "%H:%M"))) +
facet_grid(Subject ~ Affect)Nevertheless, I found some discrepancies between the data I got and the descriptions in the literature.
All analyses were conducted using the maximum available sample size: for analyses involving ESM data n=95; for all other analyses n=99.
#|label: num-of-missing-data-1
t1 <- data %>%
group_by(Subject, Affect) %>%
summarise(num_day = length(unique(Day)),
num_obser = n(),
num_missing = sum(is.na(Score)),
num_obser_valid = num_obser - num_missing) %>%
ungroup()
rmarkdown::paged_table(t1)summary(t1)
## Subject Affect num_day num_obser
## 1 : 2 Length:194 Min. :7.000 Min. :62.00
## 2 : 2 Class :character 1st Qu.:7.000 1st Qu.:66.00
## 3 : 2 Mode :character Median :7.000 Median :66.00
## 4 : 2 Mean :7.031 Mean :66.75
## 5 : 2 3rd Qu.:7.000 3rd Qu.:67.00
## 6 : 2 Max. :8.000 Max. :78.00
## (Other):182
## num_missing num_obser_valid
## Min. : 0.000 Min. :39.00
## 1st Qu.: 3.000 1st Qu.:58.00
## Median : 6.000 Median :61.00
## Mean : 6.737 Mean :60.02
## 3rd Qu.: 9.000 3rd Qu.:64.00
## Max. :29.000 Max. :73.00
##
table(t1$num_day)
##
## 7 8
## 188 6#|label: num-of-missing-data-2
t2 <- data %>%
group_by(Subject, Affect, Day) %>%
summarise(num_moment = n(),
num_missing = sum(is.na(Score)),
num_valid = num_moment - num_missing) %>%
ungroup()
rmarkdown::paged_table(t2)summary(t2)
## Subject Affect Day num_moment
## 52 : 16 Length:1364 1 :194 Min. : 5.000
## 70 : 16 Class :character 2 :194 1st Qu.: 9.000
## 87 : 16 Mode :character 3 :194 Median :10.000
## 1 : 14 4 :194 Mean : 9.494
## 2 : 14 5 :194 3rd Qu.:10.000
## 3 : 14 6 :194 Max. :20.000
## (Other):1274 (Other):200
## num_missing num_valid
## Min. : 0.0000 Min. : 0.000
## 1st Qu.: 0.0000 1st Qu.: 8.000
## Median : 1.0000 Median : 9.000
## Mean : 0.9582 Mean : 8.536
## 3rd Qu.: 1.0000 3rd Qu.:10.000
## Max. :10.0000 Max. :11.000
##
table(t2$num_moment)
##
## 5 6 7 8 9 10 11 14 20
## 10 20 138 42 188 864 98 2 2
table(t2$num_valid)
##
## 0 2 3 4 5 6 7 8 9 10 11
## 2 4 2 30 46 89 147 208 348 448 40Moment is nested into Day
\[ N_{ik(j)} = \mu +S_i + D_j + M_{k(j)} + (SD)_{ij} + \epsilon_{ik(j)} \]
where
\(N_{ijk}\): the score of negative affect for subject \(i\) at moment \(k\) in day \(j\).
\(\mu\) : the grand mean.
\(S_i\): the random effect for the subject \(i\).
\(D_j\): the random effect for the day \(j\).
\(M_{k(j)}\): the random effect of the moment \(k\) and it is nested in the specific day \(j\).
The model in Schönbrodt et al. (2022) was designed as that Day and Moment are crossed each other. However, in our data, Moment nested within Day is more reasonable because the time points of each measurements were random in each day and there are no specific meaning of those sampling points.
Follow the Schönbrodt et al. (2022) and Shrout & Lane (2012), the reliability of our model should be distinguished as between subject reliability and within subject reliability. However, these two reliabilities are defined by correlation format not by GT theory, which Schönbrodt et al. (2022) and Shrout & Lane (2012) approached (To be honest, I still cannot figure out how they built out their definition). Furthermore, in Schoenbrodt et al’s (2022) model, they assume that Moment is crossed with Day, the model has a Couple effect over subjects and item effect within each subject, and Moment effect and Item effect are fixed. So, You’ll notice a clear difference between the way I define these reliaiblities and the way they define them.
\[ \bar{N}_{i} = \frac{1}{n_D n_{M_{ij}}}\sum_{j=1}^{n_D}\sum_{k=1}^{n_{M_{ij}}} N_{ik(j)}= \mu +S_i + \frac{1}{n_D}\sum_{j=1}^{n_D}D_j + \frac{1}{n_{M_{ij}}}\sum_{k=1}^{n_ {M_{ij}}} M_{k(j)} + \frac{1}{n_D}\sum_{j=1}^{n_D}(SD)_{ij} + \frac{1}{n_D n_{M_{ij}}}\sum_{j=1}^{n_D}\sum_{k=1}^{n_{M_{ij}}} \epsilon_{ik(j)} \]
\[\begin{align} R^{bw}_{abs} &= Corr(\bar{N}_{i}, \bar{N'}_{i}) \\ &= \frac{Cov \left( \begin{matrix} \mu +S_i + \frac{1}{n_D}\sum_{j=1}^{n_D}D_j + \frac{1}{n_{M_{ij}}}\sum_{k=1}^{n_ {M_{ij}}} M_{k(j)} + \frac{1}{n_D}\sum_{j=1}^{n_D}(SD)_{ij} + \frac{1}{n_D n_{M_{ij}}}\sum_{j=1}^{n_D}\sum_{k=1}^{n_{M_{ij}}} \epsilon_{ik(j)}, \\ \mu +S_i + \frac{1}{n_D}\sum_{j'=1}^{n_D}D_{j'} + \frac{1}{n_{M_{ij'}}}\sum_{k=1}^{n_ {M_{ij'}}} M_{k(j')} + \frac{1}{n_D}\sum_{j=1}^{n_D}(SD)_{ij'} + \frac{1}{n_D n_{M_{ij'}}}\sum_{j=1}^{n_D}\sum_{k=1}^{n_{M_{ij'}}} \epsilon_{ik(j')} \end{matrix} \right) }{\sqrt{Var(\bar{N}_{i})Var(\bar{N'}_{i})}} \\ &= \frac{\sigma^2_S}{\sigma^2_S + \sigma^2_D / n_D + \color{red}{\sigma^2_M} / n_M+ \sigma^2_{SD} / n_D + \sigma^2_{\epsilon} / n_D n_M} \end{align}\]
However, the numbers of moments for each subject are not same because of the missing data. Currently, I set all \(n_{M_{ij}}\)’s equal to 10 for convenience.
From day to day or from moment to moment
\[ \bar{N}_{ij} = \frac{1}{n_{M_{ij}}}\sum_{k=1}^{n_{M_{ij}}} N_{ik(j)}= \mu +S_i + D_j + \frac{1}{n_{M_{ij}}}\sum_{k=1}^{n_{M_{ij}}} M_{k(j)} + (SD)_{ij} + \frac{1}{n_{M_{ij}}}\sum_{k=1}^{n_{M_{ij}}} \epsilon_{ik(j)} \]
\[\begin{align} R^{WS}_{interday, abs} &= Corr(\bar{N}_{ij}, \bar{N}_{ij'} | S_i = s_i) \\ &= \frac{Cov \left( \begin{matrix} \mu +s_i + D_j + \frac{1}{n_{M_{ij}}}\sum_{k=1}^{n_{M_{ij}}} M_{k(j)} + (sD)_{ij} + \frac{1}{n_{M_{ij}}}\sum_{k=1}^{n_{M_{ij}}} \epsilon_{ik(j)}, \\ \mu +s_i + D_{j'} + \frac{1}{n_{M_{ij'}}}\sum_{k=1}^{n_{M_{ij'}}} M_{k(j')} + (sD)_{ij'} + \frac{1}{n_{M_{ij'}}}\sum_{k=1}^{n_{M_{ij}}} \epsilon_{ik(j')} \end{matrix} \right) }{\sqrt{Var(\bar{N}_{ij})Var(\bar{N}_{ij'})}}\\ &= \frac{\color{red}{??????}}{\sigma^2_D + \sigma^2_M/n_M + \sigma^2_{SD} + \sigma^2_{\epsilon}/n_M} \end{align}\]
\[\begin{align} R^{WS}_{intermoment, abs} &= Corr(Y_{ik(j)}, Y_{ik'(j)} | S_i = s_i) \\ &= \frac{Cov \left( \begin{matrix} \mu +s_i + D_j + M_{k(j)} + (sD)_{ij} + \epsilon_{ik(j)}, \\ \mu +s_i + D_{j} + M_{k'(j)} + (sD)_{ij} + \epsilon_{ik'(j)} \end{matrix} \right) }{\sqrt{Var(Y_{ik(j)})Var(Y_{ik'(j)})}}\\ &= \frac{\sigma^2_D + \sigma^2_{SD}}{\sigma^2_D + \sigma^2_M + \sigma^2_{SD} + \sigma^2_{\epsilon}} \end{align}\]
Compare to Schönbrodt et al. (2022)
#|label: fitted-by-lme4
library(lme4)
neg_data <- data %>%
filter(Affect == "Neg",
Day %in% 1:7) %>%
drop_na()
neg_lme <- lmer(Score ~ (1 | Subject) + (1 | WDay) + (1 | WDay:Moment) + (1 | Subject:WDay),
data = neg_data)
summary(neg_lme)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Score ~ (1 | Subject) + (1 | WDay) + (1 | WDay:Moment) + (1 |
## Subject:WDay)
## Data: neg_data
##
## REML criterion at convergence: 43793.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.1503 -0.4552 -0.1431 0.2612 5.8739
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject:WDay (Intercept) 34.24 5.851
## Subject (Intercept) 108.01 10.393
## WDay:Moment (Intercept) 0.00 0.000
## WDay (Intercept) 0.00 0.000
## Residual 90.14 9.494
## Number of obs: 5797, groups:
## Subject:WDay, 678; Subject, 97; WDay:Moment, 80; WDay, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 15.482 1.086 14.25
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
pos_data <- data %>%
filter(Affect == "Pos",
Day %in% 1:7)
pos_lme <- lmer(Score ~ (1 | Subject) + (1 | WDay) + (1 | WDay:Moment) + (1 | Subject:WDay),
data = pos_data)
summary(pos_lme)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Score ~ (1 | Subject) + (1 | WDay) + (1 | WDay:Moment) + (1 |
## Subject:WDay)
## Data: pos_data
##
## REML criterion at convergence: 49730.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0791 -0.5708 0.0766 0.6421 3.2019
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject:WDay (Intercept) 57.993 7.615
## Subject (Intercept) 155.594 12.474
## WDay:Moment (Intercept) 3.788 1.946
## WDay (Intercept) 4.509 2.123
## Residual 262.216 16.193
## Number of obs: 5794, groups:
## Subject:WDay, 678; Subject, 97; WDay:Moment, 80; WDay, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 57.66 1.56 36.95Caution: When I ran both model estimations, lme4 displayed a warming about “boundary (singular) fit”.
| Negative Affect | Positive Affect | |
|---|---|---|
| Between subject reliability | \[ \frac{108.01}{108.01 + \frac{0}{7} + \frac{0}{10} + \frac{34.24}{7} + \frac{90.14}{7\times10}} = 0.946 \] | \[ \frac{155.59}{155.59 + \frac{4.51}{7} + \frac{3.79}{10} + \frac{57.99}{7} + \frac{262.22}{7\times10}} = 0.923 \] |
| Within subject change reliability | \[ \frac{0 + 34.24}{0 + 0 + 34.24 + 90.14} =0.275 \] | \[ \frac{4.51 + 57.99}{4.51 + 3.79 + 57.99 + 262.22} = 0.190 \] |
# Not run yet
library(brms)
neg_brm <- brm(Score ~ (1 | Subject) + (1 | WDay) + (1 | WDay:Moment) + (1 | Subject:WDay),
data = neg_data)
summary(neg_brm)
pos_brm <- brm(Score ~ (1 | Subject) + (1 | WDay) + (1 | WDay:Moment) + (1 | Subject:WDay),
data = pos_data)
summary(pos_brm)ggplot(data, aes(x = as_datetime(paste(today(), Time)))) +
geom_histogram(color = "white", binwidth = 30*60, boundary = 0) +
scale_x_datetime(name = "Time", date_breaks = "2 hours", date_labels = "%H:%M") morning: 9:00-12:00
early-afternoon: 12:00-15:00
late-afternoon: 15:00-18:00
early-evening: 18:00-21:00
late-evening: 21:00-24:00
data2 <- data %>%
filter(Day %in% 1:7) %>%
mutate(Weekend = WDay %in% c("Sat", "Sun"),
CMoment = factor(cut(as.double(Time),
breaks = c(9, 12, 15, 18, 21, 24) * 3600,
label = FALSE),
labels = c("morning",
"early-afternoon", "late-afternoon",
"early-evening", "late-evening"),
order = TRUE)) %>%
group_by(Subject, Weekend, Affect, CMoment) %>%
summarise(N = sum(!is.na(Score)),
AveScore = mean(Score, na.rm = TRUE))
rmarkdown::paged_table(data2)
data2 %>% filter(Affect =="Neg") %>%
group_by(CMoment) %>%
summarise(n_ = n(),
N = sum(N))
## # A tibble: 5 × 3
## CMoment n_ N
## <ord> <int> <int>
## 1 morning 189 853
## 2 early-afternoon 194 1375
## 3 late-afternoon 194 1515
## 4 early-evening 194 1455
## 5 late-evening 189 599\[ N^*_{ik(j)} = \mu +S_i + W_j + M^*_{k(j)} + (SW)_{ij} + \epsilon_{ik(j)} \]
where
\(N^*_{ijk}\): the new average score of negative affect for subject \(i\) at moment \(k\) in weekend \(j=1\) or not \(j=0\)
\(\mu\) : the grand mean.
\(S_i\): the random effect for the subject \(i\).
\(W_j\): the fixed effect to indicate whether the observation is in weekend or not.
\(M_{k(j)}\): the random effect of the new categorical moment $k$, including: morning, early-afternoon, late-afternoon, early-evening, and late-evening, with 5 levels. Moment* is nested in \(W_j\).
\[ \begin{align}R^{bw}_{abs} &= Corr(\bar{N}^*_{i}, \bar{N'}^*_{i} \mid W_j = w_j) \\&= \frac{Cov \left( \begin{matrix} \mu +S_i + \frac{1}{n_W}\sum_{j=1}^{n_W}w_j + \frac{1}{n_{M^*_{ij}}}\sum_{k=1}^{n_{M^*_{ij}}} M^*_{k(j)} + \frac{1}{n_W}\sum_{j=1}^{n_W}(Sw)_{ij} + \frac{1}{n_W n_{M^*_{ij}}}\sum_{j=1}^{n_D}\sum_{k=1}^{n_{M^*_{ij}}} \epsilon_{ik(j)}, \\ \mu +S_i + \frac{1}{n_W}\sum_{j=1}^{n_W}w_{j} + \frac{1}{n_{M^*_{ij}}}\sum_{k=1}^{n_{M^*_{ij}}} M^*_{k'(j)} + \frac{1}{n_W}\sum_{j=1}^{n_W}(Sw)_{ij} + \frac{1}{n_W n_{M^*_{ij'}}}\sum_{j=1}^{n_W}\sum_{k=1}^{n_{M^*_{ij'}}} \epsilon_{ik'(j)} \end{matrix} \right) }{\sqrt{Var(\bar{N}^*_{i})Var(\bar{N'}^*_{i})}} \\&= \frac{\sigma^2_S + \sigma^2_{SW}/n_{SW}}{\sigma^2_S + \color{red}{\sigma^2_M} / n_M+ \sigma^2_{SW} / n_W + \sigma^2_{\epsilon} / n_W n_M}\end{align} \]
Because Weekend is a fixed effect, I’m not sure how to build up reliability under this situation.
neg_data2 <- data2 %>% filter(Affect == "Neg")
neg_lme2 <- lmer(AveScore ~ 1 + Weekend + (1 | Subject) + (1 | Weekend:CMoment) + (1 | Subject:Weekend),
data = neg_data2)
summary(neg_lme2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: AveScore ~ 1 + Weekend + (1 | Subject) + (1 | Weekend:CMoment) +
## (1 | Subject:Weekend)
## Data: neg_data2
##
## REML criterion at convergence: 6415
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1753 -0.4418 -0.1000 0.3452 5.7003
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject:Weekend (Intercept) 14.99 3.872
## Subject (Intercept) 104.09 10.203
## Weekend:CMoment (Intercept) 0.00 0.000
## Residual 31.88 5.646
## Number of obs: 946, groups:
## Subject:Weekend, 194; Subject, 97; Weekend:CMoment, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 15.3918 1.1375 13.531
## WeekendTRUE -0.1508 0.6671 -0.226
##
## Correlation of Fixed Effects:
## (Intr)
## WeekendTRUE -0.291
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
pos_data2 <- data2 %>% filter(Affect == "Pos")
pos_lme2 <- lmer(AveScore ~ 1 + Weekend + (1 | Subject) + (1 | Weekend:CMoment) + (1 | Subject:Weekend),
data = pos_data2)
summary(pos_lme2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: AveScore ~ 1 + Weekend + (1 | Subject) + (1 | Weekend:CMoment) +
## (1 | Subject:Weekend)
## Data: pos_data2
##
## REML criterion at convergence: 7355.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0769 -0.5248 0.0424 0.5137 3.5369
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject:Weekend (Intercept) 17.325 4.162
## Subject (Intercept) 148.001 12.166
## Weekend:CMoment (Intercept) 6.134 2.477
## Residual 96.931 9.845
## Number of obs: 946, groups:
## Subject:Weekend, 194; Subject, 97; Weekend:CMoment, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 57.175 1.770 32.305
## WeekendTRUE 1.669 1.795 0.929
##
## Correlation of Fixed Effects:
## (Intr)
## WeekendTRUE -0.506| Negative Affect | Positive Affect | |
|---|---|---|
| Between subject reliability | \[ \frac{104.09 + \frac{14.99}{2}}{104.09 + \frac{14.99}{2} + \frac{0}{5} + \frac{31.88}{2\times5}} = 0.972 \] | \[ \frac{148.00 + \frac{17.33}{2}}{148.00 + \frac{17.33}{2} + \frac{6.13}{5} + \frac{96.93}{2\times5}} = 0.935 \] |
\[\begin{align} N_{it} &= \tau_{it} + \epsilon_{it} \\ &= \mu_{it} + \lambda_{T_{it}} \xi + \lambda_{S_{it}} \zeta_t + \epsilon_{it} \end{align}\]
\[\begin{align} N_{it} &= T_{i} + S_{it} \\ &= (\mu_{it} + \lambda_{T_{it}} \xi + \vartheta_i) + (\lambda_{S_{it}} \zeta_t + \epsilon_{it}) \end{align}\]
\(\mu_{i}\): the intercept (the ground mean).
\(\xi\): the common trait score.
\(\vartheta_j\): the unique trait score for \(i\).
\(\zeta_t\): the common state score.
\(\epsilon_{it} \overset{i.i.d.}{\sim} N(0, \sigma^2_\epsilon)\): the unique state score.
Common consistency: \(CCon = \frac{\lambda^2_{T_{it}}Var(\xi)}{Var(N_{it})}\)
Unique consistency: \(UCon = \frac{Var(\vartheta_i)}{Var(N_{it})}\)
\[\begin{align} N_{it} &= \tau_{it} + \epsilon_{it} \\ &= \mu_{it} + \lambda_{T_{it}}\xi_i + \lambda_{S_{jt}}O_t + \epsilon_{it} \\ O_t &= \beta_{t-1, t}O_{t-1} + \zeta_t \\ O_1 &= \zeta_1 \end{align}\]
Predictability by trait: \(Con = \frac{\lambda^2_{T_{it}Var(\xi_t)}}{Var(N_{it})}\)
Unpredictability by trait: \(Con = \frac{\beta^2_{t-1, t}\lambda^2_{S_it}Var(O_{t-1})}{Var(N_{it})}\)
I found the three kinds of variance coefficients are hard to use correlations to define them.
Level 1 (within subject)
\[\begin{pmatrix}P_{it} \\ N_{it}\end{pmatrix} = \begin{pmatrix} y_{1it} \\ y_{2it}\end{pmatrix} = \begin{pmatrix}\mu_{1i} \\ \mu_{2i}\end{pmatrix} + \begin{pmatrix}\tilde{y}_{1it} \\ \tilde{y}_{2it}\end{pmatrix} + \begin{pmatrix}\epsilon_{1it} \\ \epsilon_{2it}\end{pmatrix} \]
\[ \begin{pmatrix}\epsilon_{1it} \\ \epsilon_{2it}\end{pmatrix} \sim N_2\left(\begin{pmatrix}0 \\ 0\end{pmatrix}, \begin{pmatrix} \sigma^2_{\epsilon 11i} & \sigma^2_{\epsilon 12i} \\ \sigma^2_{\epsilon 12i} & \sigma^2_{\epsilon 22i}\end{pmatrix}\right) \]
\[ \begin{pmatrix}\tilde{y}_{1it} \\ \tilde{y}_{2it}\end{pmatrix} = \begin{pmatrix} \phi_{11i} & \phi_{12i} \\ \phi_{12i} & \phi_{22i} \end{pmatrix} \begin{pmatrix}\tilde{y}_{1it-1} \\ \tilde{y}_{2it-1}\end{pmatrix} + \begin{pmatrix}\omega_{1it} \\ \omega_{2it}\end{pmatrix} \]
\[ \begin{pmatrix}\omega_{1it} \\ \omega_{2it}\end{pmatrix} \sim N_2 \left(\begin{pmatrix}0 \\ 0\end{pmatrix}, \begin{pmatrix} \sigma^2_{\omega 11i} & \sigma^2_{\omega 12i} \\ \sigma^2_{\omega 12i} & \sigma^2_{\omega 22i}\end{pmatrix}\right) \]
Level 2 (between subject)
\[ \left(\begin{array}{c}\mu_{1 i} \\\mu_{2 i} \\\Phi_{11 \mathrm{i}} \\\Phi_{12 \mathrm{i}} \\\Phi_{21 \mathrm{i}} \\\Phi_{22 \mathrm{i}}\end{array}\right) \sim M v N\left(\left(\begin{array}{c}\gamma_{\mu 1} \\\gamma_{\mu 1} \\\gamma_{\Phi 11} \\\gamma_{\Phi 12} \\\gamma_{\Phi 21} \\\gamma_{\Phi 22}\end{array}\right), \boldsymbol{\Psi} =\begin{pmatrix} \underset{2\times2}{\boldsymbol{\Psi}_\mu} & \underset{2\times4}{\boldsymbol{\Psi}_{\mu\Phi}} \\ \underset{4\times2}{\boldsymbol{\Psi}^\top_{\mu\Phi}} & \underset{4\times4}{\boldsymbol{\Psi}_\phi}\end{pmatrix}\right) \]
\[ \left(\begin{array}{l}\sigma_{\omega_{11 i}}^{2} \\\sigma_{\omega_{12 i}}^{2} \\\sigma_{\omega_{22 i}}^{2}\end{array}\right) \sim M v N\left(\left(\begin{array}{l}\gamma_{\sigma_{\omega 11}^{2}} \\\gamma_{\sigma_{\omega 12}^{2}} \\\gamma_{\sigma_{\omega 22}^{2}}^{2}\end{array}\right), \boldsymbol{\Psi}_{\sigma_{\omega}^{2}}\right) \]
\[ \left(\begin{array}{l}\sigma_{\epsilon_{11 i}}^{2} \\\sigma_{\epsilon_{12 i}^{2 i}}^{2} \\\sigma_{\epsilon_{22 i}}^{2}\end{array}\right) \sim M v N\left(\left(\begin{array}{l}\gamma_{\sigma_{\epsilon 11}^{2}} \\\gamma_{\sigma_{\epsilon 12}^{2}}^{2} \\\gamma_{\sigma_{\epsilon 22}^{2}}^{2}\end{array}\right), \boldsymbol{\Psi}_{\boldsymbol{\sigma}_{\epsilon}^{2}}\right) \]
Total variance of data: \(Var(\mathbf{y}_{ti}) = \boldsymbol{\Psi}^2_\mu + E_i[\mathbf{T}_{\tilde{y}i}] +\boldsymbol{\gamma}_{\sigma^2_\epsilon}\)
Proof: (\(\mathbf{y}_{ti} = \begin{pmatrix} y_{1ti} \\ y_{2ti} \end{pmatrix}\) vector representation)
\[ \begin{align*}Var(\mathbf{y}_{ti}) &= Var(\boldsymbol{\mu}_i + \mathbf{\tilde{y}}_{ti} + \boldsymbol{\epsilon}_{ti} ) \\&= Var(\boldsymbol{\mu}_i) + Var(\mathbf{\tilde{y}}_{ti}) + Var(\boldsymbol{\epsilon}_{ti}) ) \\&= Part 1 + Part 2 + Part 3\end{align*} \]
Part 1:
\[\begin{align*}Var(\boldsymbol{\mu}_i)&= Var_i(E_t[\boldsymbol{\mu}_i \mid I=i]) + E_i[Var_t(\boldsymbol{\mu}_i \mid I=i)] \\&= Var_i(\boldsymbol{\mu}_i) + E_i[0] \\&= \begin{pmatrix} \psi^2_{\mu11} \psi_{\mu12} \\ \psi_{\mu12} \psi^2_{\mu22} \end{pmatrix} = \boldsymbol{\Psi}^2_\mu \end{align*}\]
Part 2:
\[ \mathbf{\tilde{y}}_{ti} = \boldsymbol{\Phi}\mathbf{\tilde{y}}_{t-1 i} + \boldsymbol{\omega}_{ti} = \boldsymbol{\Phi}^{t}\mathbf{\tilde{y}}_{0 i} + \sum_{k=1}^{t}\boldsymbol{\Phi}^{t-1}\boldsymbol{\omega}_{ti} \\Var_t(\mathbf{\tilde{y}}_{ti}) = \boldsymbol{\Phi}Var_t(\mathbf{\tilde{y}}_{t-1 i})\boldsymbol{\Phi}^\top + Var_t(\boldsymbol{\omega}_{ti}) \]
Consider \(t\) is at steady state. Let \(\mathbf{T}_{\tilde{y}i} := Var_t(\mathbf{\tilde{y}}_{ti})\), (cf., Kim & Nelson,1999, p. 27)
\[\begin{align*} \mathbf{T}_{\tilde{y}i} &= \boldsymbol{\Phi}\mathbf{T}_{\tilde{y}i}\boldsymbol{\Phi}^\top + \boldsymbol{\Sigma_{\omega i}} \\vec(\mathbf{T}_{\tilde{y}i}) &= vec(\boldsymbol{\Phi}\mathbf{T}_{\tilde{y}i}\boldsymbol{\Phi}^\top) + vec(\boldsymbol{\Sigma_{\omega i}}) \\ vec(\mathbf{T}_{\tilde{y}i}) &= (\boldsymbol{\Phi} \otimes \boldsymbol{\Phi}) vec(\mathbf{T}_{\tilde{y}i})+ vec(\boldsymbol{\Sigma_{\omega i}}) \quad (\text{by Lemma 1} \\vec(\mathbf{T}_{\tilde{y}i}) &= (\mathbf{I}-(\boldsymbol{\Phi} \otimes \boldsymbol{\Phi})^{-1} )vec(\boldsymbol{\Sigma_{\omega i}}) \\\mathbf{T}_{\tilde{y}i} &= mat((\mathbf{I}-(\boldsymbol{\Phi} \otimes \boldsymbol{\Phi})^{-1} )vec(\boldsymbol{\Sigma_{\omega i}})) \end{align*}\]
Than, we have
\[\begin{align*} Var(\mathbf{\tilde{y}}_{ti})&= Var_i(E_t[\mathbf{\tilde{y}}_{ti} \mid I=i]) + E_i[Var_t(\mathbf{\tilde{y}}_{ti}\mid I=i)] \\&= Var_i(\boldsymbol{\Phi}^{t}\mathbf{\tilde{y}}_{0 i}) + E_i[\mathbf{T}_{\tilde{y}i}] \\ &= E_i[\mathbf{T}_{\tilde{y}i}] \end{align*}\]
Part 3:
\[\begin{align*} Var(\boldsymbol{\epsilon}_{ti}))&= Var_i(E_t[\boldsymbol{\epsilon}_{ti} \mid I=i]) + E_i[Var_t(\boldsymbol{\epsilon}_{ti} \mid I=i)] \\&= Var_i(\mathbf{0}) + E_i \left[\begin{pmatrix} \sigma^2_{\epsilon11i} \sigma^2_{\epsilon12i} \\ \sigma^2_{\epsilon12i} \sigma^2_{\epsilon22i}\end{pmatrix}\right] \\ &= \begin{pmatrix} \gamma_{\sigma^2_{\epsilon11}} \gamma_{\sigma^2_{\epsilon12}} \\ \gamma_{\sigma^2_{\epsilon12}} \gamma_{\sigma^2_{\epsilon22}} \end{pmatrix} = \boldsymbol{\gamma}_{\sigma^2_\epsilon} \end{align*}\]
\[\begin{align*} R_B(\mathbf{y}) &= \frac{\Psi^2_\mu}{Var(\mathbf{y})} \\ &=\frac{Cov\begin{pmatrix} \boldsymbol{\mu_i} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{it-1} + \boldsymbol{\omega}_{it} + \boldsymbol{\epsilon}_{it}, \\ \boldsymbol{\mu_{i'}} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{i't-1} + \boldsymbol{\omega}_{i't} + \boldsymbol{\epsilon}_{i't} \end{pmatrix}}{\sqrt{Var(\boldsymbol{\mu_i} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{it-1} + \boldsymbol{\omega}_{it} + \boldsymbol{\epsilon}_{it})Var(\boldsymbol{\mu_{i'}} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{i't-1} + \boldsymbol{\omega}_{i't} + \boldsymbol{\epsilon}_{i't})}} \\ &\color{blue}{= Corr(\mathbf{y}_{it}, \mathbf{y}_{i't})} \end{align*}\]
\[\begin{align*} R_w(y_i) &= \frac{\tau^2_{\tilde{y}i}}{v(y_i)} \\ &= \frac{???}{\boldsymbol{\Phi}Var_t(\mathbf{\tilde{y}}_{t-1 i})\boldsymbol{\Phi}^\top + Var_t(\boldsymbol{\omega}_{ti}) + \sigma^2_{\epsilon_i}} \\ &=\frac{Cov\left(\begin{matrix} \boldsymbol{\mu_i} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{it-1} + \boldsymbol{\omega}_{it} + \boldsymbol{\epsilon}_{it}, \\ \boldsymbol{\mu_i} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{it'-1} + \boldsymbol{\omega}_{it'} + \boldsymbol{\epsilon}_{it'} \end{matrix} \mid i \right)}{\sqrt{Var(\boldsymbol{\mu_i} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{it-1} + \boldsymbol{\omega}_{it} + \boldsymbol{\epsilon}_{it} i)Var(\boldsymbol{\mu_i} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{it'-1} + \boldsymbol{\omega}_{it'} + \boldsymbol{\epsilon}_{it'} \mid i )}} \\ &\color{blue}{= Corr(\mathbf{y}_{it}, \mathbf{y}_{it'} \mid i)} \end{align*}\]
where \(v(y_i) = \tau_{\tilde{y}i} + \sigma^2_{\epsilon i}\) and \(\tau_{\tilde{y}i}\) is the person-specific variance for variable \(\tilde{y}^j\) and it is equal to the \(j\)th diagonal element for that variable of the covariance matrix \(\mathbf{T}_{\tilde{y}i}\) for person i.